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Abstract 

We introduce EChemDID, a model to describe electrochemical driving force in reactive molecular 
dynamics simulations. The method describes the equilibration of external electrochemical potentials 
(voltage) within metallic structures and their effect on the self consistent partial atomic charges used in 
reactive molecular dynamics. An additional variable assigned to each atom denotes the local potential 
in its vicinity and we use fictitious, but computationally convenient, dynamics to describe its equi¬ 
libration within not-simply connected metallic structures on-the-fly during the molecular dynamics 
simulation. This local electrostatic potential is used to dynamically modify the atomic electronega¬ 
tivities used to compute partial atomic changes via charge equilibration. Validation tests show that 
the method provides an accurate description of the electric fields generated by the applied voltage 
and the driving force for electrochemical reactions. We demonstrate EChemDID via simulations of 
the operation of electrochemical metallization cells. The simulations predict the switching of the de¬ 
vice between a high-resistance to a low-resistance state as a conductive metallic bridge is formed and 
resistive currents that can be compared with experimental measurements. In addition to applications 
in nanoelectronics, EChemDID could be useful to model electrochemical energy conversion devices. 

Keywords: molecular dynamics, electrochemistry, reactive potentials, electrochemical metallization 

cell 

* Corresponding author: strachan@purdue.edu 


1 



1 Introduction 


Molecular dynamics n (MD) modeling has become a key tool in many areas of science and engineering 
including chemistry, 0 0 material science 0 E] and biology [6]. Cloud computing has significantly 
simplified access to these MD simulations [ 7 ] to the point where they can now be used in undergraduate 
education [8]. While ab initio MD [9] (where atomic forces are computed from a quantum mechanical 
electronic structure calculation) plays an important role in many applications [101E] ? if s computational 
intensity limits its use to relatively small systems and short times. Therefore, dynamical simulations 
with empirical interatomic potentials represent an indispensable tool to connect electronic processes with 
the response of materials or biological systems and in determining the atomistic mechanisms that govern 
their behavior. The last decades witnessed the development of accurate interatomic potentials for a wide 
range of materials and applications including metals nain], semiconductors mus], oxides [laiTT], 
and molecular systems [HiTniiin]- More recently, reactive potentials Eunaiis] have demonstrated the 
ability to describe complex chemistry including the reaction and decomposition of high-energy density 
materials m, the description of solid/liquid interfaces [25l|26] and the mechanisms of complex catalytic 
chemical reaction m- 

A wide range of applications of technological interest require the description of electrochemical pro¬ 
cesses; examples include batteries [28], capacitors m, and electrochemical metallization (ECM) cells 
for nanoelectronics m- Despite recent progress [SUES] we lack a self-consistent approach to describe 
how externally applied electrochemical potentials affect partial atomic charges and drive electrochemi¬ 
cal reactions. In this paper we introduce the electrochemical dynamics with implicit degrees of freedom 
(EChemDID) method that describes the electrochemical driving force resulting from the application of an 
external voltage to metallic electrodes in reactive MD simulations. A significant challenge addressed by 
our method is describing the equilibration of the external potential over metallic electrodes with changing 
composition and morphology as ions electrochemically dissolve into the electrolyte or are deposited. We 
demonstrate the power of EChemDID via simulations of nanoscale ECM cells of interest for resistive 
switching non-volatile memory applications. The original version of the approach used a cluster analysis 
to determine the extent of the electrodes and enabled the first atomistic description of the operation 
of ECM cells [33]. This paper extends and validates the method, provides a detailed description of the 
underlying physics and further exemplifies its use. 

The remainder of the paper is organized as follows. Section 2 describes the EChemDID method and 
its implementation within reactive MD simulations. Section 3 verifies our implementation and validates 
the model while Section 4 exemplifies its use to describe ECM cells of interest in nanoelectronics. Einally, 
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conclusions are drawn in Section 5. 


2 The EChemDID method 

This section starts with a brief review of reactive MD simulations focusing on the self-consistent calcula¬ 
tion of partial atomic charges using charge equilibration (QEq [34]). We then introduce the EChemDID 
approach to describe the application of an electrochemical potential to metallic electrodes. 

2.1 Reactive MD and charge equilibration 

Reactive interatomic potentials like ReaxEE [21], REBO [35l|36] or COMB [23], use the concept of bond 
order to describe covalent interactions that include bond stretch, angle bending and torsional potentials 
as well as terms to penalize over and under coordination. In addition to covalent interactions, van der 
Waals interactions describe London dispersion and Pauli repulsion and electrostatics are described in 
terms of partial atomic charges. Theses charges need to be computed self-consistently based on the 
atomic structure of the system to capture charge transfer during chemical reactions. ReaxEE uses the 
charge equilibration method [34], also known as electronegativity equalization (EEM [37]), where the 
atomic charges are obtained equating the electronic chemical potential of every atoms in the system. 
The atomic electronegativity is obtained from a simple expression for the total electrostatic energy of the 
system: 


E{{qih{Ri}) = Y, (xUi + iHiqf] +YqiQjJi\Ri “ Rj\) (1) 

i ^ ' i<j 

where the second term in the right-hand side (RHS) represents the electrostatic interaction between 
ions i and j; J takes the form of shielded Coulomb interaction. The first term represents the energy 
cost associated with changing the charge of individual atoms and depends on their electronegativity (y^) 
and hardness {Hi). Given an atomic structure, the equilibrium partial atomic charges are obtained by 
equating the chemical potential of each atom: 


Xi = Xo = = ,^0 ^ Hiqi + Y QjJ i\Ri - Rj\) 


dqi 


( 2 ) 


leading to N equations with N 1 unknowns (charges for the N atoms and the chemical potential 
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Xo)- The system of equations is closed by constraining the total charge in the system: Q = qi. 

Important recent advances to charge equilibration methods include methods based on charge transfer 
between bonded atoms like the Split Charge Equilibration (SQE |38]), the addition of integer charge 
changes to describe reactions with ions dissolved in an electrolyte [31] and more recently, the extension of 
SQE and its derivation from density functional theory (DET) as proposed in the atom-condensed Kohn- 
Sham DET approximated to second order method (ACKS2 [39]). In this paper we develop EChemDID 
starting from the original formulation of QEq that remains widely used, the method is easily extensible 
to SQE. 

We are interested in simulating electrochemical reactions that originate from the application of an 
external electrochemical potential difference ($o) between two electrodes; this has the effect of changing 
the relative energy of electrons in the two electrodes by Erom the energy expression in Eq. 

we see that this can be accomplished by changing the atomic electronegativity of the atoms in one 
electrode to Xi ^ Xi ~ ^o/2 and of those in the other electrode to Xi ^ Xi + ^o/2. Changes in atomistic 
electronegativity has been used in the past to create external fields and electrostatic potentials jioisiiisg. 
The EChemDID method introduced in the next subsection, describes the process of equilibration of the 
applied electrochemical potential in non-simply-connected structure of each contact accounting for the 
loss and gain of atoms during electro-dissolution and metallization. 

2.2 EChemDID equations for electrochemical potential equilibration 

The problem we need to address now is how the external potential applied to a group of atoms in the 
electrodes propagates and equilibrates throughout the metallic structure given that their constituent 
atoms change as electrochemical reactions occur at the electrodes. We start by assigning an additional 
dynamical variable to each metallic atom that represents the local electrochemical potential in their 
vicinity: This is done in the same spirit as the local electronic temperature in two-temperature 

models designed to capture the thermal transport role of conduction electrons in metals [laES] or to 
represent the temperature of internal degrees of freedom in coarse grain simulations m. We consider that 
the electrochemical potential is externally applied to a pre-determined set of atoms located away from the 
chemically active regions, see Eigure[^(at 0 ps), and the electrochemical potential of these atoms remain 
constant throughout the simulation. We now propose a dynamical equation for the temporal evolution 
of the electrochemical potential of all atoms in the system. 

The voltage applied to a group of atoms within the metallic electrode propagates through it at the 
speed of light following Maxwell’s wave equations [45] until the entire metal is at the same electrochemical 
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potential. We are not interested in the actual dynamics of this voltage equilibration but on the much 
slower electrochemical reactions; thus we use fictitious diffusive dynamics to describe the equilibration 
process. We assume that the voltage within the metal follows: 


^ (3) 

where the dot denotes time derivative, is the Laplacian operator and k is an effective diffusivity. If 
a voltage difference is applied to electrodes that are not connected to each other, the solution of Eq. [^will 
result in the voltage being equilibrated within each electrode at a timescale that depends on the value of 
k. Furthermore, we will assume that if an atom becomes detached from the electrode its electronegativity 
should go back to its atomic value that is, the external electrochemical potential ^ should go to zero. 
The numerical solution of Eq. is performed on-the-fiy during the MD simulation using atoms as a grid 
and a local weighting function as in the eleDID method for thermal transport [43]: 


*<(() = E 






(4) 


The first term on the RHS solves the diffusion equation (Eq. |^; Rij = — Rj denotes the distance 

between atoms i and j and w{R) is a local weighting function defined below. The second term in the 
RHS relaxes the electronegativity of atoms detached from the electrodes to their atomic value. 77 is 
the relaxation rate, F{W) is a switching function that turns on the relaxation as an atom becomes 
detached from other metallic atoms and Wi is the total metallic coordination of atom i. This external 
electrochemical potential is then added to the atomic electronegativity x*(t) = Xi (0 

used for charge equilibration. 

The local weighting function takes the form: 


AT 


w{Rij) = < 



0 


if R < Rc, 
otherwise. 


(5) 


with JV a normalization constant. The range of this weighting function {Rc) represents a critical 
separation distance below which two atoms are considered as part of the same metallic cluster and, thus. 
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be at the same electrochemical potential. We define the normalization factor N from a reference lattice 
(in D dimensions, with neighbors at distance i?„) as: 


M = 


2DNat 

T.iWi 


(6) 


Nat is the total number of atoms and Wi represents the total metallic coordination of atom i in the 
reference lattice Wi = ^j^i'w{Rij). The last term in Eq. relaxes the external potential of isolated 
atoms to zero. The switching function F{Wi) increases from zero to one as the total metallic coordination 
decreases to zero as: 


F{W) 



ifW < fTo, 


otherwise. 


(7) 


Wq represents a critical metallic coordination below which an atom will evolve its electronegativity 
towards the equilibrium value. We take this value as that of the weighting function w when two atoms 
are separated by a distance close to but smaller than the cutoff Rc (Eq. [^. In practice, we take this 
value as w{0.99Rc)‘ 

Eigure exemplifies the equilibration of the chemical potential in a two electrodes setup. Each 
electrode consists of Cu atoms forming a perfect fee lattice and an applied external electrochemical 
potential of ±4 V is applied to the regions shown at time t = 0 ps. Solving the EChemDID equations 
(with atoms fixed in space in this first example) shows the equilibration of the voltage in each of the 
contacts. 

The following Sections demonstrate the coupling of EChemDID with charge equilibration that gener¬ 
ates electrostatics and with ionic dynamics to simulate electrochemical reactions. 


2.3 EChemDID and electronic transport: estimating current densities 

While the main purpose of the method is to equilibrate the electrochemical potential in metallic structures 
for electrochemical simulations, an interesting by-product of the approach is that Eqs. and i(r/ = 0) 
can be used to estimate electrical currents when a potential difference is applied across a metallic system. 
Assuming diffusive transport, a combination of Ohm’s law V4> = pj (where j is the electrical current 
and p the electrical resistivity) with the continuity equation = pVj = pq (where q is the electrical 
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1.0 ps 


10.0 ps 


Figure 1: Propagation and equilibration of the chemical potential computed with the EChemDID method. 
The arrows on the top snapshot point toward thin layers of atoms where we apply the voltage = 8 V. 
The colors represent the local electrochemical potential ^i{t) ranging from xo ~ ^o/2 (red) to xo + ^o/2 
(blue) and we use the fictitious diffusion coefficient /c = 6 A^/fs (Eq. [^. 

charge density) shows that the time derivative in the electrochemical potential in EChemDID Eq. is 
proportional to the time derivative of the free charge density as = kpq. Let’s consider a continuous 
region of a metallic wire of cross sectional area A sandwiched between the regions and the 
time derivative of the total charge Qo in Q. can be expressed as: 





Qn 


— f ^ ^ 

ien 






( 8 ) 


with Qi an atomic volume. Erom Eq. (assuming r] = 0) Qq can be written as the sum of two terms 
representing the current flowing to/from and 


Qn = - 
P 




|i?, 


^3 \ 


E 


E. 

iG^^ 3 


I ^ij I 


w{Rij)fti 


(9) 


Einally, under steady state, the total current density from the EChemDID calculation can be obtained 
as: 




iG^^ 3 — 1 




wiRij) = E E 






w{Rij) 


( 10 ) 
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The equality of the currents coming in and out of the region of interest offers a simple verification of 
the numerical implementation. 

2.4 EChemDID implementation 

The EChemDID method has been implemented as an external LAMMPS [46] user-package m including a 
LAMMPS “compute” that solves the Laplacian in Eq. [^and a “fix” that integrates the voltage diffusion in 
time (Eq. [^. The whole implementation is consistent with the parallel scheme employed in LAMMPS. 
The modifications are independent to the actual LAMMPS program, only minor modifications to the 
reax/c package [48] allow switching from regular QEq (with constant electronegativity) to dynamical 
electronegativity in EChemDID. In the following MD simulations, the interactions between atoms are 
described by ReaxEE for Cu-Si02 from the Ref. [33]. While EChemDID has been fully implemented 
following the model presented in the previous section, a simplified version will be applied in the following 
ignoring the relaxation term and we will refer to the numerical solution of the Laplacian as Eq. @(^ = 0). 
This last approximation assumes that the particles retain memory of their previous electronegativity. 

3 Implementation verification and EChemDID validation 

We first verify our implementation of the EChemDID equilibration using atoms as a grid by analyzing 
the ID propagation of the electrochemical potential through a metallic electrode. We then validate the 
approach by computing the electric fields generated by the atomic charges induced by the difference in 
potential between two electrodes. 

3.1 Electrochemical potential equilibration and electrostatics 

We start by verifying the EChemDID equilibration of voltage by comparing the numerical solution with 
an analytical one in a ID transport over a contact of length L with one end set to +4 V. The initial 
voltage through the sample is zero and the analytical solution for the diffusion equation leads to voltage 
profiles as a function of time given by: 


with Xn = ny/kir/L. 


;r. +C>0 ^ 

4>o 1 

— > —e smX 

^ n L 


TT 


( 11 ) 
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Figure 2: Propagation of voltage along a copper contact computed with EChemDID (Eq. compared 
to the analytical solution (Eq. [n]). The chemical potential profile is represented every 10 fs from 0 to 
100 h {k = 6 A^/fs and $o/2=+4 V). 

Eigure compares the numerical solution (red circles) with the analytical one (lines) at different 
times; the EChemDID numerical integration uses a timestep of 0.05 fs. The excellent agreement verifies 
our numerical solution of the equation Eq. Eor this example, the Cu atoms are fixed in space; coupling 
Eq. I^with molecular dynamics will be demonstrated in Section]^ 

Having demonstrated that EChemDID equilibrates the electrochemical potential within metallic elec¬ 
trodes we now show how changing the atomic electronegativity affects the QEq-derived atomic charges 
and, thus, creates electric fields. We do this first for a simple parallel plates capacitor setup separated 
by vacuum. The application of an electrochemical potential 4>o between infinite, parallel electrodes sep¬ 
arated by vacuum of thickness d should generate a constant electric field between them equal to 4>o/d 
and a surface charge density equal to a = ^o/deo (with eo the relative permittivity). After solving the 
EChemDID equations and performing charge equilibration we compute the magnitude of the electric field 
between the electrodes using a probe charge that scans space. The resulting force F on the probe atom of 
charge q is proportional to the electric field as F = qE. Eigure shows 2D maps of the resulting electric 
field generated by two types of electrodes. The electric field created by the fiat electrodes represented in 

is spatially constant whereas the roughness induced by the triangular electrode showed in|^ results in 
electric field enhancement at the tip of the electrode. In both cases depicted in Eigure the amplitude 
of the electric field agrees with the applied voltage with approximately 5% of error. 

3.2 Driving force for electrochemistry 

In order to quantify the driving force for electrochemistry we compare the energy of a series of atomic 
configurations under different voltages. We move a probe consisting of a Cu atom through a Inm-thick 
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Figure 3: Maps of the electric field between two electrodes under a potential of 8 V: a) the electric 
field generated by two perfectly flat electrodes (separated by a distance d = 1.37 nm, resulting in the 
voltage $0 ~ 7.6 V); b) the electric field generated by a triangular shaped electrode versus a perfectly 
flat electrode (separated by a tip-to-fiat distance d = 0.98 nm, resulting in the tip-voltage ~ 8-3 V) 
highlighting the effect of field enhancement at the tip of the patterned electrode. 



Figure 4: Energy profile of a system composed by a copper probe atom (located in z) dissolved in 
Si02 sandwiched between two copper electrodes at different voltages. The actual voltage applied to the 
electrodes is shown in the legend and the effective voltage is shown in parenthesis. 


layer of crystalline silicon dioxide sandwiched between two copper electrodes. At each location of the 
probe atom, we relax the structure via energy minimization and plot in Figure the resulting energy 
landscape. The voltage applied between the parallel metallic electrodes acts on the ion dissolved in the 
dielectric and the average slope of the energy landscape is given by: dEjdz — g4>o/d where d—X nm is 
the separation between electrodes and g=0.35 e the (averaged) charge of the dissolved ion. Therefore, we 
can extract an effective voltage from the total energies as a function of ion position obtained from the 
simulations. Figure [^compares the extracted effective (in parenthesis) to the applied voltages. We find 
a reasonable agreement between the two considering the large fluctuations originating from the rugged 
underlying atomistic energy landscape. 
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4 Simulating electrochemical reactions using EChemDID 


We now focus on full EChemDID simulations where Eq. (77 = 0) is solved together with a reactive MD 
simulation. To demonstrate the approach we simulate the operation of nanoscale ECM cells of interest in 
nanoelectronics [49]. These resistance-switching devices consist of an electro-active electrode (typically 
Cu or Ag) and an electro-inactive one (often W or Pt) separated by a solid electrolyte (amorphous Si 02 
in our case). When a voltage of the appropriate polarity is applied, the atoms in the active electrode 
dissolve into the electrolyte and are field-driven towards the inactive electrode. Eventually, this process 
leads to the formation of a metallic conductive filament between the two electrodes, with the subsequent 
decrease in cell resistance. 

4.1 ECM simulation details 

The simulation cell is composed by two metallic electrodes separated by amorphous Si 02 generated 
following an annealing procedure described in Ref. EO] and similar to the setup presented in Ref. [33] . 
The active electrode is made of Cu and so is the inactive one, which is described as a rigid body to 
capture its electrochemical inertness. In order to mimic the roughness of real interfaces we randomly 
deleted 50 % of the Cu atoms part of the first layer at the interfaces. The structure is then relaxed via 
energy minimization and equilibrated using isobaric/isothermal MD for 50 ps. The ECM cell is 4x4 nm^ 
in cross section where we apply periodic boundary conditions and the separation between electrodes is 
approximately 1.0 nm. 

To simulate the operation of the cell, we impose a constant potential ±4>o/2 to a thin layer of atoms 
(^ 0.2 nm) at each boundary of the electrodes (in the direction perpendicular to the interfaces) resulting 
in a voltage of 8 V, characteristic during the forming phase of this type of device m- The MD equations 
of motion are solved with a timestep of 0.5 fs within the canonical ensemble at 300 K. The voltage 
diffusion equation (Eq. is integrated 10 times per MD timestep (for numerical stability) and the 
diffusion constant is chosen to be equal to 4 A^/fs with a cutoff distance Rc =0.4 nm (Eq. [^. Eor the 
purpose of demonstration, we use EChemDID with we the coefficient set to 0; a movie of the ECM 
simulation is included in the Supplementary Material. 

4.2 Resistance Switching 

Eigurej^ shows atomistic snapshots at various times during the switching of the ECM cell, for clarity only 
the copper atoms are shown; the bottom panel shows the time evolution of the corresponding switching 
state of the cell. The color on the top row of snapshots denotes the local electronegativity (blue is xo-4 
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V and red Xo+4 V) and the bottom row uses color to denote partial atomic charges. We see that QEq 
correctly localizes charges at the interfaces between the copper electrodes and the amorphous silica with 
positive ions on the active electrode and negative at the inactive one. The low electronegativity in the 
active electrode leads to positive partial charges which facilitates its dissolution into the dielectric; this 
can be formally thought of as an oxidation reaction: Cu ^ +(5e“. Similarly, the electronegativity of 

the dissolved copper ions increases as they become closer to the inactive electrode (within the cutoff Rc), 
lowering their partial charge until they bond to the electrode and become metallic following chemical 
reduction, formally: Cu‘^+ + Se~ Cu. Throughout the simulation, EChemDID maintains constant 
electrochemical potential within each electrode even as their morphology changes dramatically during 
the operation of the device. As shown above, the charged electrodes generate an electric field throughout 
the dielectric that acts as the driving force for the diffusion of copper ions inside the solid electrolyte 
from the active toward the inactive electrode during the forming or programing phase. 

EChemDID also provides insightful information regarding the formation of a conducting filament. 
To study this process we use two approaches: i) a distance-based cluster analysis to determine whether 
the two electrodes are connected by a Cu filament as was done in Ref. [33] and ii) using the current 
computed from EChemDID Eq. The distance-based cluster analysis uses a cutoff of 0.4 nm and 
we assign a value of 1 to the signal when the two electrodes are bridged (low resistivity state: LR), 0 
otherwise (high resistivity state: HR); the red line in Eigure (bottom scale) shows a running average 
(with a time-window of 25 ps) of the resulting signal for clarity. Therefore, a signal between LR and HR 
corresponds to switching states where the cell oscillates between ON and OEE within the time-window 
of the average. The current represented on the bottom scale of Eigurej^ (blue) has also been averaged 
with a time-window of 5 ps. 

Erom the cluster analysis, we see that a filament was established for a short duration just before 0.5 
ns and a more stable one forms at time 1.4 ns; this second filament remains for the duration of the run. 
The atomic configuration of the first bridge consists of a linear single-atom chain (inset b) and bridges 
the two electrodes for approximately 0.1 ns. As in our prior work [33] single-atom chains are metastable 
and do not survive long timescales. The second, more stable, filament (inset c) is slightly thicker with 
a zig-zag shape. We show in the insets of Eigurej^ (b and c) the time-averaged coordination number of 
the atoms in the filament. While the atoms in the first bridge exhibit an averaged coordination number 
between 2 and 3, typical of a linear configuration, the atoms in the zig-zag chain have higher coordination 
number suggesting a stronger connection. 


The current computed via Eq. 10 correlates with the distance-based cluster analysis however ad- 
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a) t = 0 ns b) t = 0.5 ns c) t = 1.4 ns 



0 0 0 5 10 1.5 

time InsJ 


Figure 5: Snapshots of the ECM cell at various times during switching and corresponding switching state 
of the device computed from Eq. (blue, in arbitrary units) or by a distance-based cluster analysis 
(red, HR/LR = high/low resistance state). The colors on the top panels represent the electronegativity 
on each atom ranging from 0.8 V (blue) to 8.8 V (red). The bottom panels show partial atomic charges 
ranging from -0.3e (red) to +0.3e (blue). The a-Si02 between the electrodes has been hidden for clarity. 
We show on the insets in b and c zoomed-in representations of the filament and the local (time-averaged) 
coordination number of the constituent atoms. The coordination numbers in b have been time-averaged 
from 0.4 to 0.5 ns and the numbers in c have been averaged between 1.4 to 1.6 ns. 


ditional informations can be extracted from its amplitude. Contrary to the binary state of the cluster 
analysis, the amplitude of the current corroborates the strength of the connection, hence the quantity of 
current that can flow through the filament. We see on Eigurej^ a peak in the current approximately 3 
times larger than in b in agreement with the inset c showing a zig-zag configuration of the bridge. We 
demonstrate a direct connection between the atomic structure of the filament and the current leaving and 
entering the electrodes computed from the dynamical equilibration of the EChemDID method. This last 
analysis is the qualitative in silico equivalent to atomically-controlled quantum conductance experiments 
proposed recently [5^ . 
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5 Discussion and conclusions 


In summary, we introduced a novel approach to describe electrochemical reaction using reactive molecular 
dynamics simulations. EChemDID describes the equilibration of electrochemical potential throughout 
non-simply connected metallic structures and this potential shifts the energetics used to obtain atomic 
charges in charge equilibration. Thus, the application of an external potential leads to atomic charges 
that can drive electrochemical process. We demonstrated that the driving force for electrochemistry is 
accurately described by EChemDID. In addition, the electrochemical potential equilibration is performed 
on-the-fly during the MD simulation and thus adjusts to changes of composition and topology of the 
electrodes. As atoms dissolve into the electrolyte their electronegativity evolves back to its original value 
or keeps memory of previous voltage. The implementation in LAMMPS is particularly straightforward 
and the method can be easily extended to more sophisticated charge equilibration methods. 

EChemDID is not without limitations. The ions dissolved in the electrolyte could trap electrons and 
change their charge state; such processes are currently ignored in EChemDID. Recently, Miiser and Dapp 
proposed a method to handle integer charge states and used it to explore redox reactions in a model 
battery; such an approach could be incorporated into EChemDID. In addition, while EChemDID enables 
an estimation of electronic currents, they do not affect atomic dynamics; processes like Joule heating and 
electromigration are ignored in this first version of the model. 

The demonstration of the use of EChemDID to simulate the resistance switching in ECM cells shows 
the power of the approach which is generally applicable with reactive and non-reactive force fields (as 
long as dynamical charges are computed with QEq) and can be used to study electrochemical processes 
with atomistic detail. There is growing interest in ab initio descriptions of electrochemistry [53], and 
EChemDID can help bridge such methods toward full device simulations. 
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